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Abstract 

We propose a stochastic lattice gas model to describe the dynamics of two 
animal species population, one being a predator and the other a prey. This 
model comprehends the mechanisms of the Lotka-Volterra model. Our anal- 
ysis was performed by using a dynamical mean-field approximation and com- 
puter simulations. Our results show that the system exhibits an oscillatory 
behavior of the population densities of prey and predators. For the sets of 
parameters used in our computer simulations, these oscillations occur at a 
local level. Mean-field results predict synchronized collective oscillations. 
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I. INTRODUCTION 



The dynamics of interacting species has deserved renewed interest since a simple math- 
ematical model had been proposed independently by Lotka and Volterra ||. The 
Lotka- Volterra model predicts an oscillatory temporal evolution of a predator-prey system 
and has been widely used as a jumping-off place for models of many organism societies 
||-[5|]. However, some of its appealingly simple assumptions fail in taking into account rele- 
vant features of such systems. For instance, an individual is generally affected by its local 
environment and not by the global density of each specie. The importance of space, con- 
sidered as a limiting resource which would enhance competitive coexistence, has also been 
pointed out 0. It is in this context that lattice gas versions of the model caught recently 
new attention PHID] as a suitable way to address these questions and simulate practical 
examples. 

We propose a model consisting of a system of two types of interacting particles residing 
in the sites of a lattice. One type of particle represents a prey and the other a predator. 
Each site can be either empty (0), occupied by one prey (X) or occupied by one predator 
(Y). The system evolves in time according to a stochastic irreversible dynamics with local 
interactions similar to those occurring in contact process models [|3~T| [X3T and also in lattice 
gas models describing chemical reactions fn|,[To]|. The local interactions considered are: 



predators can be spontaneously annihilated, Y — > 0, prey can be autocatalytically created, 
+ nX — > (n + 1)X, and predators can also be autocatalytically created at the expense of 
prey, X + nY —>■ (n + 1)Y. In this stochastic treatment of the system, with local evolution 
rules, fluctuations inherent to a many-body irreversible interacting system are taken into 
account from a microscopic point of view. 

We have analyzed the states of the system as a function of the creation rates of prey and 
predators and of the death rate of predators considering square lattices and one dimensional 
lattices. Our computer simulation and mean-field analysis indicate that the system can 
exhibit four kinds of states: two absorbing states, one being a prey-absorbing state and the 
other a vacuum-absorbing state; and two active states, one presenting an oscillatory behavior 
in the population densities and the other characterized by constant stationary population 
densities. 

The layout of the paper is as follows. We start in Sec. |TJ defining the model and setting 
up its master equation. In Sec. |T| we use mean field approximations to solve the system 
of equations associated to the master equation. There is a qualitative description of the 
different type of solutions that we have found for the one-site and the pair approximations. 
Sec. |V| describes how computer simulations were performed and the results that we obtained. 
A summary and concluding remarks are stated in Sec. [V]. 

Finally, let us make a short digression on terminology. In this work we use many times 
the terms 'prey' and 'predator'. Eventhough inspired in mathematical ecology, they are only 
generic names. Actually, words like 'hosts' or 'healthy cells' and 'parasites' or 'infected cells' 
may be more accurate to mimic a biological process. Or we may simply call them particles 
of type A and B. 
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II. THE MODEL 



Consider a lattice of N sites which can be either empty, occupied by a prey or occupied 
by a predator. At each time step a site is randomly chosen. For that site, we denote by 
n a( n b) the number of nearest neighbors occupied by prey (predators), and by ( the total 
number of nearest neighbors. Within this model three cases need to be considered: 

(a) The site is empty and it becomes either occupied with a prey with a probability 
Pa n a/( or else it remains empty with probability 1 — p a n a /(- 

(b) The site is occupied by a prey, which is either replaced by a predator with a probability 
Pb n b/( or eke remains as a prey with probability 1 — Pb n b/C- 

(c) The site is occupied by a predator, which is either vacated with a probability p c or 
else remains as a predator with probability 1 — p c . 

The markovian process defined above involves three parameters: p a , pb , andp c , which are 
associated to three subprocesses. Subprocess (a) describes the birth of prey, subprocess (b) 
the death of prey and simultaneous birth of predators, and subprocess (c) the spontaneous 
death of predators. 

The state of the system is represented by a = (<7i, a"2, <Jn) where <Ji = 0, 1 or 2 accord- 
ing to whether the site % is empty, occupied by a prey or occupied by a predator. 

Let P(a, t) be the probability of state a at time t and let Wi{a)/r be the probability per 
unit time of a cyclic permutation of variable <7j. That is, if <7j = 0, 1 or 2, then wi is the 
transition probability to cr, = 1, 2 or 0, respectively. The evolution of P(u, t) is governed by 
the master equation, 

d N 

r-P(a, t) = t) - Wi (a)P(a, t)}, (1) 

at - =1 

where the state denoted by a 1 is obtained from state a by an anticyclic permutation of the 
variable <Tj. 

According to the local rules of the model, defined above, we have: 



Wi(a) =p a -^5(<7j +(5 , 1) if Oi = 0, (2a) 
s <5 
1 

Wi(a) = p b - ^25{a i+s , 2) if o { = 1, (2b) 
s 5 

Wi(a) = p c if Oi = 2, (2c) 



where the summation is over the ( nearest neighbor sites and 5(x, y) is the Kronecker delta. 

By rescaling the time the process is found to be invariant under the transformation 
p a —> ap a , pb — > apb, and p c — > ap c where a is a positive constant. Hence we restrict 
ourselves to a set of parameters such that they satisfy the condition p a + pb + p c = 1 • This 
is automatically satisfied by writing p a = 1/2 — p — c/2, p b = 1/2 + p — c/2, and p c = c. All 
the figures have been made in terms of parameters p and c. 

The average of the state function f(a) is defined by: 

</(*)>=E/(*) p (M)- (3) 
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From Eqs. (0) and the master equation (P, the time evolution of < f(a) >is 
d N 

A < M >= E < IfM ~ MW?) >, (4) 

where the state denoted by V is obtained from a by a cyclic permutation of the variable av 
Let us make the following definitions: 

Pi(a) = < 5(<Ti,a) >, (5a) 
Pijiptp) = < 5{<Ti, a)8(aj,P) >, (5b) 
Pijk{a(3-i) = < 6(ai,a)5((Tj, (3)5(a k ,j) >, (5c) 



where a, /3, and 7 can take any one of the values 0, 1, or 2 . Using Eqs. (|2]) we may write 
the time evolution equations for some of these probabilities as 

rjU(l) =p a W P iM 01 ) -p&^E^W 12 ), (6a) 

T jP l (2)= Pb Wp l M^) -p c mi (6b) 

T ± Pij (01) = - Pa l £ (P i+ ^.(10l)+^(0l))+^ E (PaPi j>j+ s(0l0) - p b Pi j>j+ 5(012)) 

+p c ^(21), (6c) 
7-^(12) = i E (Pa^+^(102)-p 6 Pi +5 ,y(212))+p^ E ^W*( 112 ) 

-(P6T+Pc)^i(12), (6d) 

r^(02)=p c i^(22)-p a i ^ Pi+6^102) + - E p 6 P^(012) -p c J^(02). (6e) 

These expressions constitute a hierarchic system of equations. The time evolution of the 
one-site correlations Pi(a) involve the two-site correlations Pij(aj3), the time evolution of the 
two-site correlations involve the three-site correlations Pijk(af3'j), and so on. A truncation 
of this hierarchic system, which we consider in the next section, is the starting point for 
obtaining approximate solutions to the problem. 



III. TRUNCATION APPROXIMATION 



In order to obtain approximate solutions of Eqs. (y) we use a truncation scheme fl6}-|19 
The simplest truncation scheme is obtained by writing the probability of a cluster of sites 
as the product of the probability of each site. A truncation of higher order, that is of order 
n > 1, consists in writing any correlation in terms of correlations of order n and less than 
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n. Consider a cluster of m > n sites and denote it by C. Let A and B be the sets of points 
in the 'core' and 'boundary' of cluster C, respectively. The 'core' A is chosen to have n — 1 
sites. The conditional probability V(B\A) is approximated by the product Ilies ^(^l-^)- 
Therefore the probability V{C) of cluster C is given by 

V{C) = V{A)V{B\A) « P(A) J] P(<|A) = P(A) J] (7) 

i&b i&b 1 y A ) 

where V(i, A) is the probability of the cluster of n sites formed by site i and the sites of A. 



A. One-site approximation (n = 1) 

In this approximation one obtains a closed set of equations for the one-site correlation 
Pi(a). This is accomplished by writing any two-site correlations Pij(a/3) as the product 
Pi(a)Pj(/3). Three one-site correlations are present, Pi(0), -Pi(l), and P&(2), but only two 
can be independent. We choose Pi(l) and Pi(2) so that Pi(0) = 1 — Pi(l) — P(2). The time 
evolution of these variables are then 

rjU(l) = pJpPm E Pi+si 1 ) - Pb^m) E p i+s( 2 )> (8a) 
Tj t P l (2)=p b ^P t (2)Y^P l+ 5(l)-PcP(2). (8b) 

Since we seek for homogeneous solutions we write P»(l) = x, Pi (2) = y and Pj(0) = 
1 — x — y. Eqs. (g) become 

T—x = p a (l-x-y)x-p b xy 1 (9a) 
= Pfc^Z/ - Pcy- (9b) 

at 

These equations admit two trivial fixed points, (x,y) = (1,0) and (x,y) = (0,0), which 
correspond to the prey-absorbing and vacuum-absorbing states, respectively. The linear 
stability analysis reveals that the latter is always unstable and the former is stable in region 
I of the phase diagram (fig. [I]). Below the critical line, defined by c = (2p + l)/3, this 
point becomes unstable giving rise to a nontrivial stable fixed point for which < x < 1 and 
< y < 1, namely an active state. Two regions can be distinguished according to the way the 
nontrivial fixed point is approached. In region II, the fixed point is an asymptotically stable 
node (real eigenvalues) whereas in region III it is an asymptotically stable focus (complex 
eigenvalues). This implies the emergence of damped oscillations in the population densities 
of the system. The line that separates regions II and III is given by p a p c = ^Pb{pb — Pc), 
where p a , p b and p c are defined in terms of p and c in Sec. |I[ 



B. Pair approximation (n = 2) 

This approximation consists in writing the three-site correlations in terms of two-site 
and one-site correlations. In this case equation (|7|) leads to the following expression for the 
probability of a cluster of three sites, 
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Wr) - ^»), (10) 

where sites z and are nearest neighbors of site j. 

We also seek for spatially homogeneous and isotropic solutions of Eqs. @. In this case 
we may drop the indexes in Pi(a) and Pjj(a/3). We then have three one-site correlations 
P(a), a = 0,1,2 and nine two-site correlations P(a/3), a,/3 = 0, 1,2 . However, only five 
of them are independent. We choose them to be -P(l) = x, the prey density, P(2) = the 
predator density, P(01) = u, P(12) = i>, and P(02) = w. The equations for these variables 
are 

d 

r—x = p a u-p b v, (11a) 
d 

r-j^y = Pbv - p c y, (lib) 

T ~dt U= 7 \ P<1 ~ Pb— ) + PcV - Pa-U, (lie) 




T ~dt V = 7 — \ Pa ~ + Pb ~ \~PbT v ~Pc v i ( lld ) 

T d\ w = ~ — ~ ■ ~ p "~ ) " ~ ir] ■ (iif') 

where z = P(0) = 1— x — y, q = P(00) = z — u — w, r = P(ll) = x — u — v , and 
s = P(22) =y-v-w. 

We have obtained analytically the stationary solutions of the system of Eqs. (|TTD. Due 
to their somewhat lengthy expressions, we thought more convenient to make a qualitative 
discussion here. Their explicit dependence on the parameters p and c for a general ( is given 
in Appendix A. 

For the case ( = 2, the system admits only one stationary solution which is the prey- 
absorbing state. 

We have analyzed these equations for ( = A. They have two trivial fixed points: 
(x,y, z,u,v,w) = (1,0,0,0,0) and (x,y, z,u,v,w) = (0,0,0,0,0) that correspond to the 
prey-absorbing and vacuum-absorbing states respectively. The linear stability analysis re- 
veals that the former becomes unstable in a line of critical points which corresponds to the 
curve separating regions I and II in fig. £| The region where the prey-absorbing state is 
stable is denoted by I. Below this critical line, the system of equations displays a non trivial 
fixed point, namely a stationary solution, with < x < 1 and < y < 1. Fig. || shows 
the densities of prey, predators and empty sites as a function of c for the case p = O.As c 
decreases three types of fixed points are obtained, which correspond to regions II, III, and IV 
of fig. 0. Regions II and III differ in the way the stable fixed point is approached. Region II 
has an asymptotically stable node whereas region III shows an asymptotically stable focus. 
At the transition line, fixed points are asymptotically one-tangent stable focuses. This is 
due to the emergence of complex eigenvalues associated to the dominant eigenvector and 
implies the emergence of damped oscillations in the population densities of the system. 

For sufficiently small values of c a Hopf bifurcation takes place at a critical line. Inside 
region IV limit cycles are present as shown in fig. |j. These results can also be seen in fig. |5|, 
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where the real component 7 and the squared imaginary component u 2 of the dominant 
eigenvalue associated to the focus are plotted against c for p = . We also observe the 
behavior uj ~ ^J~c for sufficiently small values of c . 



IV. SIMULATION 



Numerical simulations were performed in square lattices with periodic boundary condi- 
tions. Unless another size is specified, all the figures in this section correspond to simulations 
in lattices of 100 x 100 sites. 

Each run started with an initial configuration of prey, predators and empty sites placed 
randomly in the lattice. 

Simulation results are summarized in the phase diagram of fig. |6|. 

Region I represents prey-absorbing states. Increasing the predator survival probability 
(decreasing c), the system undergoes a second order kinetic phase transition towards active 
stationary states that are represented by region II in the phase diagram, see fig. [7|. These 
states are characterized by short transients of less than 1000 MCS. An active state in the 
phase diagram means that the lifetime of the system is much longer than 30000 MCS. The 
states in region II have a constant mean value of the population densities and our analysis 
does not indicate any oscillatory behavior in this region. This behavior can be also seen from 
the graph of the spectral density S(w) of their temporal evolution, shown in fig. |8|. Spectral 
densities were obtained using the Fourier transform (FT) of temporal samples j^0[ : 



where f(n) denotes the stationary population density at time n. To avoid a divergence in 
S(w) at w = we subtracted to each temporal sample its mean value, thus obtaining a 
zero average process. Spectral densities in fig. |8| and fig. || were smoothed by an averaging 
process over nearest neighbors. 

As c is decreased for a fixed value of p, a transition takes place from region II to region 
III, consisting in the emergence of temporal oscillations in the population densities of the 
system. The oscillatory behavior is corroborated by the clear maximum in the spectral 
density, as seen in fig. || The degree of oscillation is measured by the eight of the maximum 
value of the spectral density S(w mauX ). We define a transition point as that in which the 
low modes of the spectral density (lim^^o S(w)) become as relevant (have the same value) 
as its maximum S(w m!ix ). The transition line between regions II and III was determined 
following this criterion. Increasing c in region II we get closer to the second order kinetic 
phase transition (KPT). As a consequence, long range fluctuations become highly relevant 
and predominance of the low modes is to be expected. For this reason, care should be taken 
about the reliability of the right end of this line (where regions I, II and III meet in fig. [)]). 
Since the line is so close to the second order KPT, low modes become important before 
^(wmax) suffers a significant decrease, leading to an overlap in S. 




(12) 



Wj = 




J = l,N 



N 
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Simulation results indicate that oscillations become more evident as c decreases, i.e. the 
maximum of S(w ) increases as c — > 0. Also w max decreases in this limit, for p = and 
c< 1 : Wmax ~ with /3 < 1. This behavior is in good qualitative agreement with the pair 
approximation predictions. 

Increasing c, w max also grows, reaching a limit value until the oscillation disappears. For 
p = this value corresponds to a period of ~ 250 MCS. 

Region IV in fig. |B| represents vacuum-absorbing states. Using lattices up to 480 x 480 
the vacuum-absorbing region reduces its height with the lattice size of the system. This is in 
agreement with the pair-approximation predictions that the vacuum-absorbing state occurs 
only when c = 0. But it still remains to be studied the behavior of this region as L — > oo, 
in particular for the case p > (where p a < Pb). 

In regions I to IV, the system evolves towards an active or absorbing state independently 
of the initial densities of prey and predators. This is not the case for the parameter set of 
the grey region in fig. In this region, depending on the initial densities, the system will 
reach an oscillating state or will become trapped in either of its two absorbing states. In the 
last case, the time it takes for the system to become trapped is much less than the transient 
of any simulation resulting in an active state. A possible way to treat both this and the 
finite size problem (previous paragraph) is mentioned in Sec. M. 

In order to analyze, from the simulations, whether the system presents oscillatory syn- 
chronized collective states or if the oscillations occur at a local level, we studied the standard 
deviation of the temporal samples 

a = ^< (p- <p>y>. (13) 

For sets of values of the parameters p and c inside the region III, we have verified that 
increasing the lattice size L of the system a decays as 1/ \^N, where N = L 2 is the number 
of sites in the lattice. However it still remains to be seen what happens for values of c << 1 . 
In the hypothetical case that a tends to a finite value as N — > oo, the system should display 
synchronized collective oscillatory states which should correspond to the region IV in fig. |2|, 
as found in the pair approximation analysis . 

V. CONCLUSIONS 

We have studied the dynamics of two competing species from a stochastic point of view 
developed by using a lattice gas model. A regular square lattice was considered where 
each site can be empty, occupied by one prey or one predator. The system evolves in time 
according to irreversible stochastic local rules. The system presents two absorbing states: a 
prey-absorbing state and a vacuum-absorbing state. We have found, by means of a mean 
field approximation and by computer simulations, that the system also presents active states 
where prey and predators are continuously being created and annihilated. These states can 
be of two types: a stationary one and another where populations of the two species oscillate 
in time. Mean field results present two kinds of oscillating states: one in which oscillations 
are damped and another consisting of a stable limit cycle of these quantities. The transition 
from one to the other takes place through a Hopf bifurcation. 

Computer simulations also exhibit a prey-absorbing state, an active state with stationary 
populations and another active state which exhibits oscillations in the population densities. 
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These results are in qualitative agreement with mean field predictions. In the two regimes 
where the probability of spontaneous death of a predator and prey reproduction are very 
small, simulations are difficult because the system eventually becomes trapped into either 
the prey-absorbing or the vacuum-absorbing state, even for large systems. This problem 
could be overcome by using another procedure in which absorbing states are avoided mod- 
ifying the model in order to include a small probability of spontaneous prey creation (of 
order 1/N). Another alternative could be to protect the very last individual of each specie 
against extinction, avoiding completely the two absorbing states. We are still testing these 
procedures. 

In the region where simulations are conclusive we have observed that oscillations occur at 
a local level revealing that the oscillatory behaviour displayed by real prey-predator systems 
might be a local phenomenon. By this we mean that there is a finite correlation length A. 
However we do not discard the possibility of having synchronized collective oscillations in 
the system (A — > oo) in the limit where the predator death probability is very small. This 
speculation is based on the fact that mean-field results indicate that a limit cycle is stable 
in this regime. This analysis will be object of a future work. 
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APPENDIX A: 



As mentioned in Sec. |III B| we will write the stationary solutions of the Eqs. (pT|). 



Defining the auxiliary functions: 




u = —(B-VBi-4AC), 
A=(p-l)(p+l), 

B = (7 + a) 09 + 1) + (/3 - l)(a + a) + 2/3/3, 
C = (7 + a) (a + a) + 2f3a, 
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the solutions are: 



a 2u . i . . . 

x = (l + -u + - a — — Al 

v c a + a-(f3+l)u' v ' 

y = -xu, (A2) 

c 



(A3) 
(A4) 

w = {(3u — a)(l — x — y). (A5) 



f = "^^5 (A4) 



10 



REFERENCES 

* e-mail: javier@uspif.if.usp.br 
** e-mail: tcastro@uspif.if.usp.br 
[1] A. J. Lotka, Proc. Nat. Acad. Sci. U.S.A. 6, 410 (1920). 

[2] V. Volterra, Legon sur la Theorie Mathematique de la Lutte pour le Vie (Gauthier- 
Villars, Paris, 1931). 

[3] N. S. Goel, S. C. Maitra, and E. W. Montroll, Nonlinear Models of Interacting Popula- 
tions (Reviews of Modern Physics Monographs, 1971). 

[4] G. Nicolis and I. Prigogine, Self- Organization in Nonequilibrium systems (Wiley- 
Inter science, New York, 1977). 

[5] H. Haken, Synergetics (Springer- Verlag, New York, 1983). 

[6] R. Sole, J. Bascompte, and J. Vails, J. Theor. Biol. 159, 469 (1992). 

[7] K. Tainaka, and Y. Itoh, Europhys. Lett. 15, 399 (1991). 

[8] K. Tainaka, J. Phys. Soc. Jpn. 57 , 2588 (1988). 

[9] K. Tainaka, Phys. Rev. Lett. 63, 2688 (1989). 
[10] H. Matsuda, N. Ogita, A. Sasaki, and K. Sato, Prog. Theor. Phys. 6, 1035 (1992). 
[11] T. E. Harris, Ann. Probab. 2, 969 (1974). 

[12] T. M. Ligget, Interacting Particle Systems (Springer- Verlag, New York, 1985). 
[13] R. Dickman, Phys. Rev. B 40, 7005 (1989). 

[14] R. M. Ziff, E. Gulari, and Y. Barshad, Phys. Rev. Lett. 56, 2553 (1986). 

[15] T. Aukrust, D. A. Browne, and I. Webman, Phys. Rev. A 41, 5294 (1990). 

[16] H. Mamada, and F. Takano, J. Phys. Soc. Jpn. 25, 675 (1968). 

[17] R. Kikuchi, Phys. Rev. 81, 988 (1951). 

[18] R. Dickman, Phys. Rev. A 34, 4246 (1986). 

[19] T.Tome, and M. J. de Oliveira, Phys. Rev. A 40, 6643 (1989). 

[20] N. G. Van Campen, Stochastic Processes in Physics and Chemistry (North-Holland, 
Amsterdam, 1981). 



11 



FIGURES 



FIG. 1. Phase diagram in c-p plane according to the one-site approximation. Region I repre- 
sents prey-absorbing states. In region II active states are asymptotically stable nodes and in region 
III active states are asymptotically stable focuses. 

FIG. 2. Phase diagram in c-p plane according to the pair approximation. Region I represents 
prey-absorbing states. In region II active states are asymptotically stable nodes and in region III 
active states are asymptotically stable focuses. Region IV corresponds to limit cycle solutions. 

FIG. 3. Population densities of prey, predators and empty sites according to the pair approxi- 
mation as c varies for p = 0. 

FIG. 4. A typical limit cycle. In this case p = and c = 0.017. Curves were plotted by 
iterating the pair approximation eqs. for two different initial conditions A and B. Both trajectories 
are counterclockwise oriented. 

FIG. 5. Real component 7 and squared imaginary component u) 2 of the dominant eigenvalue 
associated to the focus of the pair approximation eqs. , p = while c varies through regions II, III 
and IV. 

FIG. 6. Phase diagram in c-p plane obtained from simulations for a square lattice of L = 100. 
Region I represents prey-absorbing states, region II non-oscillant active states and region III oscil- 
lant active states. Region IV consists of vacuum-absorbing states. In the grey region no reliable 
information can be drawn due to the strong initial configuration dependence. 

FIG. 7. Population densities of prey, predators and empty sites according to simulations as c 
varies for p = 0. 

FIG. 8. (a): Temporal evolution of the population densities of prey and predators for p = 
and c = 0.1795 (region II). (b): Spectral density of the prey population density (similar graphs are 
obtained for predators and empty sites). 

FIG. 9. (a): Temporal evolution of the population densities of prey and predators for p = 
and c = 0.022 (region III), (b): Spectral density of the prey population density (similar graphs are 
obtained for predators and empty sites). 
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